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We study the three-fermion problem within a ID model of a Feshbach resonance in order to gain 
insight into how the Fulde-Ferrell-Larkin-Ovchinnikov(FFLO)-like state at small negative scattering 
lengths evolves into a Bose-Fermi mixture at small positive scattering lengths. The FFLO state 
possesses an oscillating superfluid correlation function, while in a Bose-Fermi mixture correlations 
are monotonic. We find that this behavior is already present at the three-body level. We present 
an exact study of the three-body problem, and gain extra insights by considering world lines of a 
path-integral Monte Carlo calculation. 
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I. INTRODUCTION 

Trapped ultracold clouds of fermions such as ^Li pro- 
vide unique insights into the superfluidity of neutral 
fermions and have opened up new directions for inquiry. 
By considering the three-body problem, here we theoreti- 
cally address the properties of a one-dimensional (ID) su- 
perfluid gas of spin-imbalanced fermions (where > n^) 
when the interactions are tuned via a Feshbach reso- 
nance. We flnd a change in symmetry of the ground- 
state wave function as a function of system parameters, 
and connect this symmetry change with properties of the 
many-body state. Our conclusions come from (i) the 
scattering lengths calculated from an exact solution of 
the 3-body problem and (ii) the off-diagonal elements 
of the pair density matrix calculated with path-integral 
Monte Carlo. In the latter formulation the symmetry 
change in the wave function emerges from a competition 
between two classes of topologically distinct imaginary- 
time world lines. Our conclusions are relevant to experi- 
ments on ^Li atoms trapped in an array of very elongated 
traps, formed from a two-dimensional optical lattice [1]. 
When such a lattice is sufficiently strong, one has an ar- 
ray of independent ID systems, and experiments probe 
ensemble- averaged quantities including the momentum 
distribution of pairs. 

Similar experiments in three dimensions (3D) have 
demonstrated a crossover between BCS superffuidity of 
loosely bound pairs to a Bose-Einstein condensation 
(BEC) of molecules, finding particularly rich physics 
(mostly involving phase separation) when the gas is spin 
polarized [2] . One dimension brings a new set of phenom- 
ena, driven by quantum fiuctuations and the topology of 
the Fermi surface. 

Of particular interest, Fermi surface nesting in ID 
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stabilizes [s*] a version of the Fulde-Ferrell-Larkin- 
Ovchinnikov (FFLO) phase in the spin-imbalanced gas 
[4]. FFLO phases, which occupy an extremely small re- 
gion of the 3D phase diagram [5], are characterized by 
a coexistence of magnetic and superfluid order, typically 
coupled together with a spin-density wave. An intuitive 
example is given by a quasi- ID spin-imbalanced BCS su- 
perfluid, where one flnds an array of 7r-domain walls in 
the superfluid order parameter, with the excess unpaired 
atoms residing near the nodes [6]. At higher polariza- 
tions the domain walls merge, and the order parameter 
becomes sinusoidal. We are interested in the truly ID 
limit, where there is no long range order: Instead, one 
can introduce an operator b{x) which annihilates a pair 
at position flnding the analogy of FFLO state is that 
(6t(x)6(0)) - cos(27rnFx)/|x|^ where uf = Ti^ — n^ is the 
density of excess fermions and the exponent S depends on 
interactions [7]. 

When the interactions are weak, a sufficiently dilute 
and cold gas of ^Li atoms in an elongated trap (with 
transverse dimension d = ^/K/mIJ±) can be modeled as 
a ID Fermi gas interacting through a short-range ID po- 
tential [8] . This mapping requires that the 3D scattering 
length is negative with \a\/d <^ 1, and both the thermal 
energy ksT and the chemical potential /i are small com- 
pared to the transverse confinement energy huj±. Like 
Refs. [Qi, ilO|], we will consider stronger interactions. The 
breakdown of the mapping onto a ID Fermi gas is illus- 
trated by the situation where the 3D scattering length 
is small and positive, hence producing a deeply bound 
molecular state. The correct description of the unpolar- 
ized system in this limit is a weakly interacting gas of 
these bosons: a model which is not equivalent to a ID 
gas of fermions with point interactions. 

If one spin imbalances the system in this BEC limit, 
one does not produce a FFLO state, but rather the excess 
fermions only mildly perturb the bosonic pairs, and the 
correlation function (^b\x)b{0)) ^ l/\x\^ is monotonic 
prj . Here we study the three-body problem to address 
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FIG. 1: (Color online) Cartoon depictions of the physics of 
Eq. {[1]) in the BEC (left) and BCS (right) hmits. (Top) Sym- 
metry of Bose wave function: in the BCS limit the wave func- 
tion changes sign whenever a pair passes a (spin-up) fermion. 
(Middle) Depiction of lattice model which is used for devel- 
oping intuition about Eq. (1). (Bottom) Typical world lines 
illustrating interaction of a boson (heavy line) and fermion 
(thin line) with space along the horizontal axis and imaginary- 
time along the vertical axis. 



the key question of how a spin-imbalanced gas evolves 
between this fluctuating "BEC" limit and the fluctuating 
"BCS" limit already described. How does the correlation 
function go from monotonic to oscillatory? We find that 
in the three-body problem the transition occurs due to a 
level crossing. 

To this end, we consider the minimal ID model of a 
Feshbach resonance [fH, [l3| , which can capture the rele- 
vant physics, 



k,a k 



4m 
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where L is the length of the system and c|.^, Ck^a{b\, 
hk) are fermionic(bosonic) creation/annihilation opera- 
tors. The parameter g describes the coupling strength 
between the bosonic and fermionic channel and v is the 
detuning with u ^ oo {u ^ — oo) being the BCS (BEC) 
limit. We will use units in which h'^/m = 1. [23| 



II. QUALITATIVE STRUCTURE 




FIG. 2: (Color online) The dimensionless ID scatter- 
ing lengths hg/a — cis/aO'^^^ for the symmetric (solid red 
line) /antisymmetric (dashed blue line) channel plotted vs the 
dimensionless detuning z> = vjg^^^ . The dotted (dashed 
dotted) line is the asymptotic result for as, as — 3Wo^ 
(as = (3/2)z//^2) in the BCS limit (BEC limit); (cf [T^). 
(Inset) Sum as + aa (solid line) crosses zero at i> —0.635, 
marking the change in symmetry of the ground state, (c)-(e) 
Lowest-energy symmetric (solid line) /antisymmetric (dashed 
line) wave function fs/aix) = L~"^/^ J^^g e^^^ fs/a,Q^ in a hard- 
wall box of size L ^ 160/^^^^, where x represents the rel- 
ative separation of the boson and fermion. Left to right: 
i> = —1,-0.635,1. (f)-(h) Wave function near the origin. 
Finite range of the effective interaction is apparent from 
the nonsinusoidal shape of / for small x. (i)-(k) Reduced 
density-matrix p{x,x') defined in the text before Eq. (|6|) for 
P — 100/^^/^ caculated with QMC. Blue/red represents pos- 
itive/negative weight. Quadrants with predominant posi- 
tive/negative weight are labeled with "-^"/"-'\ 



Figure [T] shows a cartoon depiction of the lattice ver- 
sion of this model. One can represent the model in terms 
of two ID channels, represented as the legs of a ladder. 
Fermions move on the lower leg, while bosons move on 
the upper. As shown at (1) and (2), pairs of fermions can 
hop from the lower leg to the upper, becoming a boson 
and vice versa. 

In the BCS limit, u ^ the atoms mainly sit on the 
lower leg, making virtual transitions to the bosonic leg. 
These virtual transitions lead to a weak local attraction 



between fermions, U = —g^jv. The figure on the bottom 
right illustrates typical world lines for three fermions. 

In the BEC limit, —v^ the atoms mainly sit on 
the upper leg. They make virtual transitions to the lower 
leg. As illustrated at (3), a boson cannot make a virtual 
transition if an excess fermion sits at that location. This 
leads to a repulsive interaction between the bosons and 
fermions of strength g^ jv. Unlike the BCS limit, the 
world-lines of the fermions and bosons cross. 
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III. WAVE FUNCTIONS 



IV. QUANTUM MONTE CARLO (QMC) 



To gain insight into how this symmetry change occurs, 
we study the eigenstates of Eq. ([1]) for the case of three 
particles. Mora et al. [lo| carried out a similar study 
for a more sophisticated model of fermions confined to a 
harmonic waveguide. The simpler nature of our model, 
which only includes the most relevant degrees of freedom, 
makes the physics more transparent. 

We study what the symmetry of the ground state is 
as a function of the dimensionless parameter z> = v j g^l'^ . 
Given that the three-body wave function can be written 
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we ask what the symmetry of fx is under switching the 
relative position of the boson and the fermion (i. e. K ^ 
—K). We find that the ground state / switches from 
odd (consistent with FFLO) to even (consistent with a 
Bose- Fermi mixture) as v is increased from large negative 
values. 

To arrive at this result, we integrate out the three- 
fermion part of the wave function [14], deriving an inte- 
gral equation for the two-particle wave function /x. 



£(Q,E)/q = -^^ 



K' 



K''^ ^ QK' ^ Q'^ - E' 



(3) 



where 

C{Q, E) = 3QV4 ^i^-E- ^V(2V3QV4-^) (4) 

For details, see Appendix [Al 

The low-energy symmetric and antisymmetric scatter- 
ing states have the form V^s(x) ex sin[/c(|x| — as)] and 
tlJa{x) (X sin[A:(x + sign(x)aa)] for large \x\. By impos- 
ing hard- wall boundary conditions, f{x = ±i^/2) = 0, 
one sees that the ground state will be symmetric when 
as > —aa^ and antisymmetric otherwise. Figure [2fa) 
shows these scattering lengths as a function of reveal- 
ing that the symmetry of the wave function changes at 
z> ^ —0.635 , where the two solutions are degenerate. 
Figures [2fb)-(h) shows the structure of the lowest en- 
ergy symmetric and antisymmetric wave functions with 
these boundary conditions. Note that on the BCS side 
of resonance, where — > ^s, the Bose- Fermi interac- 
tion cannot be described W a local potential, rather it 
is a more general kernel [10]. The off-diagonal nature 
of the interaction allows the system to violate the stan- 
dard theorem that the ground-state wave function of a 
nondegenerate system has no nodes. The level crossing 
between the states of differing symmetry suggests one 
of several scenarios for the many-body system, with the 
most likely candidates being a first-order phase transition 
or a crossover. Similar behavior was seen by Kestner and 
Duan [m in their investigation of the three-body prob- 
lem in a 3D harmonic trap. 



We developed a QMC algorithm to calculate thermo- 
dynamic quantities in this model and to give alternative 
ways of thinking about the underlying physics. We cal- 
culate the thermal density matrix 

p(x, x') = Z-^ Tr \e-^^b\x)c\{0)c^{0)b{x') \ , (5) 



where b{x) = Y.k^'^^^^k, c^(^) = 
/3 = l/Zc^T, and Z is the partition function. Figures 
[2ji)-(k) shows a density plot of this correlation function. 
The FFLO phase is distinguished from the Bose-Fermi 
mixture by the sign of p in the top left and bottom right 
quadrants. The boundary between these behaviors oc- 
curs roughly where — = ^s- 

Considering first the fermionic sector, with two spin- 
up and one spin-down fermions, we discretize imaginary 
time into Af slices, writing 



/It 2t i . It 2t i. o\ _ 

p{x^^XjQ-^Xj^] Xq , Xq ,Xq; P) — 2^ 



^dxfdx^e-'^ 



(6) 



as integrals over the positions of the up-spins x^^ and 
the down-spin at imaginary times tj = jf3/JV, with 
discretized action S. For appropriately chosen this 
expression converges to the exact thermal expectation 
value as A/" oo. Two separate boundary conditions 
account for the fermionic statistics: Jj has 
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and j\f — 



while has x^ — 



It 
.It 



2t J 2, _ ^ 

X and j\f — *^ 
The integrals are performed by a Monte Carlo algorithm, 
treating as a probability measure. Details of our 
choice of discretized action and the resulting Monte Carlo 
rules are given in Appendix [Bl 

While path-integral QMC techniques are well estab- 
lished [16], the present situation is novel because two 
fermions can bind and form a boson. We implement 
this feature by introducing extra variables that record 
the slices at which two fermions are bound, and requir- 
ing that when two fermions are bound (say x]^ and xj) 
then their positions must be equal. The moves in our 
Markov process are as follows: moving a particle in one 
time slice, binding two unbound fermions of opposite spin 
into a boson, and unbinding two fermions. In all cases 
the probabilities of the move in slice j only depends on 
the positions at time slices j — 1 and j + 1. Sampling 
new positions from a Gaussian centered about weighted 
average of the particle's position in the previous and last 
slice optimizes the acceptance rate. As described in Ap- 
pendix [Bl we find the rules summarized in Table |I] and 
illustrated in Figs. [3fa)-(d), which let ([6]) converge to 
the exact density matrix as A/" ^ oo. Specifying these 
Markov rules is equivalent to specifying S. 
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TABLE I: Gaussian sampling widths and Metropolis acceptance rule, A = min(l, e ^^Tr/Tf), for moves in Figs. [3] (a)- (d). 
Moves for bead sampled from a Gaussian of width ap centered about Xj] while the reverse moves Xj Xj sample 

a Gaussian of width gr. 

Move Sampling midpoint Xj Sampling width cff Sampling width cfr e~^^TR/TF 

(a) Fermion ^.i+i'^.i-i y^A^ V^A^T^ 1 

(b) Boson "^'+^;"^-^ 1 

(c) Close^open V^A^ V^AT/^ exp(zy At) / AtVStt At 



(d) Zip^unzip 2^.-1 ySfTl VaTTI exp ( !/Ar + '•'^"+^3^7+' ) /V2 




FIG. 3: Illustrative moves in our QMC algorithm. Fermions 
are designated by thin lines with arrows representing the spin, 
bosons by thick lines, and moving beads are white: (a) Mov- 
ing a fermion, (b) moving a boson, (c) opening/closing, and 
(d) zipping/unzipping, (e) Crossing of same-spin fermions is 
always canceled by an equal weight path of opposite sign, (f) 
Bosons enable paths with both negative and positive weight 
that do not cancel. 



Since the density matrix involves adding up terms 
with different signs, at low temperatures or large parti- 
cle numbers the efficiency can suffer; this is the "fermion 
sign problem." For three particles the variance remains 
small enough that we can produce reasonably accurate 
results with the algorithm already described. To make 
further improvements, we use of the fact that paths can- 
cel when world lines for identical fermions cross in ID, a 
well-known technique for eliminating the sign problem in 
ID [2^]. For example. Fig. ^e) illustrates two paths for 
which has the same value, but which contribute to p 
with opposite signs. We therefore throw away both sets 
of paths. In a purely ID system of fermions one could 
thereby eliminate all paths with one sign or the other, 
depending on the relative ordering of the particles at the 
beginning and end. Here the cancellation is incomplete. 
Figure [3l^f) illustrates paths of opposite sign which have 
no term of the opposite sign to cancel. When the ex- 
changes are dominated by paths with positive weights 
(such as the RHS of Fig. 3(f)) one has a Bose- Fermi 
mixture; otherwise one has an FFLO-like state. 



V. REALIZATION/DETECTION 



We studied the simplest model for the BEC-BCS 
crossover of spin-polarized fermions in harmonic waveg- 
uides, a many-body system realizable by spin imbalanc- 
ing the array of ID tubes created in [1]. In such an ex- 
periment one could distinguish FFLO from a Bose-Fermi 
mixture by either using an interferometric probe [l^ 
or measuring the pair momentum distribution, e.g. by 
sweeping to the BEG side followed by time-of-flight ex- 
pansion. The signature of the FFLO phase is a peak 
at finite momentum q = ttuf set by the density of ex- 
cess fermions uf = — This peak should 
be absent in a Bose-Fermi mixture with monotonically 
decaying superfluid correlations. Another probe, based 
on correlations in the atomic shot noise after time-of- 
flight expansion, has been suggested in [19]. Additionally, 
there has recently been effort in studying the BEC-BCS 
crossover in few-body clusters [^O]. By creating ensem- 
bles of elongated clusters one can directly realize and 
study the three-body system considered here: tuning in- 
teractions using a photoassociation or a Feshbach reso- 
nance [13, EH. 
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Appendix A: Solution of the 3-body problem 

Applying the Hamiltonian Eq. [T]to |^) gives a pair of 
coupled Schrodinger equations 



9 if, 



k-\-K/2 



f-k+K/2) 



2VT 



^K"^ ly - E j /k - 2-^ ^K/2+fc,3K/4-fc/2 







0. 



In the first equation we antisymmetrized the second term 
in k to ensure manifest antisymmetry of gK,k' Eliminat- 
ing gx^k from the last equation gives an integral equation 
for fx 

37^2,,, n 9'^ \^ fK- fk-K/2 ^ ^ , ^ 

4^ +^-^J^^= L^FTS^VI^ ^^^^ 



performing 



the 



integral 



After 

JZ,dk/[2Tr{k^ + 3KyA-E)] = 1/(2^^3X^4 -E) 
this simplifies to (we relabeled K to Q) 



C{Q,E)fQ = -^J2 



K' 



L ^ K'^^ QK' ^Q^-E 



(A2) 



w\t\iC{Q,E) = 3QV4+i/-^-^/V(2v/3QV4 - E). The 
solution Eb of the equation C{0,Eb) = is the two-body 
bound state energy [l2|,[25| and C{Q, Eq) = is solved 
by Eq = 30 V4 + Eb. Equat ion ()A2[) can be converted 
into a Lippmann-Schwinger equation for the scattering 
amplitude F{Q^K) using the ansatz [10], 

/q = 27r5{K - 0) + iF{Q, K) ^ 1/{K ± Q + ze), (A3) 



which gives 

2iF{Q,K)K = -V{Q,K, K) 

(A4) 

where we introduced the effective boson fermion potential 

V{Q,K',K) = ^^^^ Ek){K'^ + QK' + Q^-Ek) 

At low energies (small momentum il^[26|) one has 
F{Q,K) -1 + iKas + iQaa ^ where as{a 
the scattering length for the symmetric (antisymmetric) 
channel. To extract the low-energy scattering proper- 
ties, we numerically solve the integral equation Eq. (|A4[) 
at fixed small K and compute the scattering length from 
the limits 

as = lim lm[F(K, K) + F(-K, K)]/2K, (A6) 
aa = ]\m^lm[F{K,K) - F{-K,K)]/2K. (A7) 



Appendix B: Derivation of the path integral action 
and Monte Carlo rules 

The partition function Z corresponding to the Hamil- 
tonian [Eq. ([1])] can be expressed as a path integral. The 
path-integral formulation is useful both as a computa- 
tional tool, but also provides insights from a different 
point of view. We formulate the path integral in real 
space (position basis) and imaginary time, < r < /3. 
We discretize imaginary time into discrete steps Ar = 
P/JV^ where JV is the Trotter number. The path in- 
tegral is equivalent to Eq. ([T]) in the limit A/" ^ 00, 
which is taken by extrapolating our numerical results to 
Ar 0, following the high- accuracy method of Schmidt 
and Lee [22| . 

To construct the path integral, we start from the imag- 
inary time propagators for free fermions, 

Gf{xa,x'a';Ar) = {0\c,{x)e-^^ ^ cl{x')\0) 

= exp (- '^^XO / 

(Bl) 

and for free bosons, 

Gb{x,x';I^T) = {0\b{x)e-^^"b\x')\0) 

( \x-x'\^ ^ \ I ^— 
= exp ( - ' ' J/ At) 

(B2) 

These propagators are illustrated in Fig. |4j To represent 
the interaction, we weigh a vertex [Fig.HJ^d)] with —g Ar. 
The action S for a path configuration, used in Eq.( [6j), 
is then given by the negative of the log of the product 
of the propagators and interaction vertices that make up 
the path. 

To sample the path, we use the Metropolis algorithm, 
in which the acceptance of a move is given by 



^ = min(l,e-^^Ti^/Ti.), 



(B3) 



where Tp and Tr are the forward and reverse probabili- 
ties of attempting a particular move. For example, con- 
sider the move illustrated in Fig.[3l^c), where a bead on a 
bosonic path is split to form a fermionic pair (a bubble). 
For the forward move, we sample the two fermion posi- 



tions and from a Gaussian of width crp 



(a) 



(b) 



A 



Gf{x],x' ]-/\t) Gf{x Ix' i-/\r) Gb{x,x'-/\T) -g Ar 

FIG. 4: Propagators and interaction vertex for disretized 
path-integral model of Eq. ([T]). (a) Spin- up free fermion, 
Eq. (|B1|) . (b) spin-down free fermion, Eq. (|B1|) . (c) free bo- 
son, Eq. (|B2|) . and (d) interaction vertex with extra weight 
-^Ar. 
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centered about Xj = (xj+i + Xj_i)/2 where Xj-i and sample the recombined boson position Xj from a Gaus- 
Xj^i are the stationary boson positions immediately be- sian of width gr^J^tJI centered about x^ We find 
fore and after the sampled slice. For the reverse move, we 



Tr _ G/(x^-+i,xJ)G/(xJ,x^_i)G/(x^-+i,xJ)G/(4,^i-i)(^Ar)2 e '^-r (27ra|)~^ 

CC J I ^ I cc ^. cc J- I ^ 

g e ^^^(27rcr|.)-i (B4) 



exp(i^ Ar) 



72 AtVStt Ar 



This rule and the rules for the other moves illustrated in Figs. [3fa)-(d) are summarized in Table HI 
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